#
# PlotSimpleResult.R
#
# Replication file for
# Hill, Seth J. and Gregory A. Huber. 2018. ``On The Meaning of Survey Reports of Roll Call `Votes','' American Journal of Political Science.
#

rm(list=ls())
if (.Platform$OS == "windows") {
  # Set working directory in location of script.
  .doit <- function() { # only works with R.exe; trap errors if using Rscript.exe
  frame_files <- lapply(sys.frames(), function(x) x$ofile);frame_files <- Filter(Negate(is.null), frame_files)
  PATH <- dirname(frame_files[[length(frame_files)]]); setwd(PATH) ; rm(PATH,frame_files)
  }
  try(.doit(),silent=T)
}

library(foreign)
library(pscl)
library(data.table)
library(RColorBrewer)
palette(brewer.pal(n=9,name="Set1")[c(1:5,7:9)]) 
library(car)
par(cex.main=1.2)

#
# Function to create slopeplot and votes with dems figures.
#
.doFigures <- function(use.weights,use.screened,first.data.file,first.fig.name,second.data.file,second.fig.name) {
  # Function to create slopeplot and votes with dems figures.
  #
  # %%%%%%%%%%%%%%%%%%%%%%%%
  # Data outsheeted from AnalyzeDemogsAndCreateRecodes.do
  # %%%%%%%%%%%%%%%%%%%%%%%%
  #
  DT <- fread(sprintf("Tables/%s",first.data.file))

  # Call in bill names.
  source('03_99_billNames.R')
  bill.names <- .billNames()
  # Create billname to merge to DT.
  bill.names[,billname := gsub("\\.","",real.short)]
  # Only care about unique match of billname to short.
  DT <- merge(DT,unique(bill.names[,c("billname","short"),with=F]),by="billname",all.x=T)
  # Merge in Senate vote splits.
  bill.votes <- .billVotes()
  DT <- merge(DT,bill.votes,by="short",all.x=T)

  #
  # %%%%%%%%%%%%%%%%%%%%%%%%
  # Comparison of information treatment to standard
  # support and Senate vote by bill and party.
  # %%%%%%%%%%%%%%%%%%%%%%%%
  #

  .slope <- function(x,y1,y2,pch=19,lab="D",col=2) {
    # Create a sloping line connecting two points.
    y1 <- unlist(y1)
    y2 <- unlist(y2)
    x0 <- x - 0.5
    x1 <- x + 0.5
    segments(x0=x0,x1=x1,y0=y1,y1=y2,lwd=2)
    points(x=c(x0,x1),y=c(y1,y2),pch=pch,col=col,cex=1.5)
    text(x1,y2,pos=4,labels=lab)
  }

  .slopeSupportPlot <- function(DT,lvar="voteStandard",rvar="voteInfo",tvar=NULL,labs=c("D","R"),member.votes=T) {
    # Make slope plot connecting support in two conditions.
    # Arguments.
    #  DT -- data.table of data.
    #  lvar -- variable measuring left hand support.
    #  rvar -- variable measuring right hand support.
    #  tvar -- variable with t ratio to display at y=-0.1.
    #  labs -- labels for two parties.
    #  member.votes -- plot MC vote support?
    par.old <- par(mar=c(4.1,3.1,3.6,1.1))
    palette(brewer.pal(n=6,name="Set1"))

    # Open plot.
    plot(x=c(1,DT[,length(unique(short))]*3),y=c(0,1),type='n',ann=F,axes=F);axis(2,las=2)
    title(ylab="Proportion supporting bill",line=2.1)
    abline(h=0.5,col='gray')
    # Plot for each bill.
    x <- 1
    for (bl in DT[,sort(unique(short))]) {
      # Dem slope.
      .slope(x=x,y1=DT[short == bl & pid3 == "Dem",lvar,with=F], y2=DT[short == bl & pid3 == "Dem",rvar,with=F], lab=labs[1], col=2)
      # Rep slope.
      .slope(x=x,y1=DT[short == bl & pid3 == "Rep",lvar,with=F], y2=DT[short == bl & pid3 == "Rep",rvar,with=F], lab=labs[2], col=1, pch=15)

      if (member.votes & rvar != "votec_split") {
        # Dem member votes.
        points(x=x+1.5,y=DT[short == bl,pdem.yea][1],col=2,pch=1,cex=1.5)
        # Rep member votes.
        points(x=x+1.5,y=DT[short == bl,prep.yea][1],col=1,pch=0,cex=1.5)
      } else if (member.votes & rvar == "votec_split") {
        # Chamber votes.
        points(x=x+1.5,y=DT[short == bl,chamber.yea][1],col='black',pch=18,cex=1.5)
      }
      if (!is.null(tvar)) {
        text(x=x+.5,y=0.01,labels=sprintf("t=%1.1f",DT[short == bl,tvar,with=F][1,]),pos=1)
      }
      # Label x-axis.
      axis(ifelse(x%%2==0,1,3),at=x+.5,labels=paste(strwrap(bl,15),collapse="\n"),tick=F,cex.axis=.8,line=ifelse(x%%2==0,2,-1))
      abline(v=x+2,lty=2)
      x <- x+3
    }
    par(par.old)
  }

  pdf(sprintf("figures/%s",first.fig.name),width=10,height=6)
  .slopeSupportPlot(DT,tvar="did_tInfo")
  dev.off()

  # Calculate average splits in each condition.
  splits <- unique(DT[,c("billname","pdem.yea","prep.yea"),with=F])
  splits[,abs.diff := abs(pdem.yea - prep.yea)]
  cat("Average vote split among Senators:",splits[,mean(abs.diff)],"\n")
  splits <- unique(DT[pid3 %in% c("Dem","Rep"),c("billname","voteInfo","voteStandard","pid3"),with=F])
  splits[pid3 == "Dem",voteInfo := -1*voteInfo]
  splits[pid3 == "Dem",voteStandard := -1*voteStandard]
  diffs <- splits[,lapply(.SD,sum),.SDcols=c("voteInfo","voteStandard"),by="billname"]
  cat("Average vote split b/w Ds and Rs, standard condition:",diffs[,mean(abs(voteStandard))],"\n")
  cat("Average vote split b/w Ds and Rs, info condition:",diffs[,mean(abs(voteInfo))],"\n")

  #
  # Re-call data.
  # 
  cat("\nCalling in individual data for other figures...\n")
  DT <- fread(sprintf("RawData/%s",second.data.file))
  setnames(DT,gsub("v1","V1",names(DT)))
  # Join recodes.
  dmg <- data.table(read.dta("RawData/Huber_and_Hill_Roll_Call_Experiment-recodes.dta"))
  DT <- merge(DT,dmg,by="V1",all.x=T,all.y=F)
  rm(dmg)
  if (use.weights) {
    # Join in post-stratification weights.
    wts <- data.table(read.dta("RawData/RakedPewWeights.dta"))
    DT <- merge(DT,wts,by="V1",all.x=T,all.y=F)
    DT[is.na(weight), weight := 1] # make Sen weights 1
    rm(wts)
  } else {
    DT[,weight := 1]
  }
  cat("Diff in diff models run excluding pure independents.\n")
  DT[pid3_DtoR != "Oth/Ind", dem := 1*(pid3_DtoR == "Dem")]
  DT[,table(pid3_DtoR,dem,exclude=NULL)]
  DT[assignedCondition != "Senator",info := 1*(assignedCondition=="Info")]
  DT[,table(assignedCondition,info,exclude=NULL)]


  #
  # %%%%%%%%%%%%%%%%%%%%%%%%
  # Distribution of vote-with-Dems-Senators
  # rate by treatment assignment.
  # %%%%%%%%%%%%%%%%%%%%%%%%
  #

  # Determine which position on each bill is "with Dems"
  # by whether Senate Dem yea rate >= Senate Rep
  # yea rate.
  bill.votes[,dem.yea := 1*(pdem.yea >= prep.yea)]
  # Grab real.short name.
  bill.votes <- merge(bill.votes,unique(bill.names[,c("short","real.short"),with=F]),by="short")
  bill.votes[,real.short := sprintf("bill.%s",real.short)]

  # For each individual, determine on each bill if they
  # voted "with Dems." 
  for (bl in bill.votes[,real.short]) {
    DT[,sprintf("vwdem.%s",bl) := 1*(DT[[bl]] == bill.votes[real.short == bl,dem.yea])]
  }

  # Aggregate number votes "with Dems" for each individual.
  DT[,tot.votes := rowSums(!is.na(DT[,grep("vwdem",names(DT)),with=F]))]
  DT[,votes.with.Dems := rowSums(DT[,grep("vwdem",names(DT)),with=F],na.rm=T)]
  DT[,table(votes.with.Dems,tot.votes)]
  DT[tot.votes > 4,prop.with.Dems := votes.with.Dems/tot.votes]

  # Collapse votes with Dems to 2-bins.
  DT[,votes.with.Dems2 := recode(votes.with.Dems,"0:1='0-1'; 2:3='2-3'; 4:5='4-5'; 6:7;'6-7'; 8:9='8-9'; 10:11='10-11';12='12'")]
  DT[,votes.with.Dems2 := factor(votes.with.Dems2,levels=c("0-1","2-3","4-5","6-7","8-9","10-11","12"))]
  #DT[,votes.with.Dems2 := factor(floor(votes.with.Dems/2),levels=1:7,labels=c("0-1","2-3","4-5","6-7","8-9","10-11","12"))]
  DT[,table(votes.with.Dems,votes.with.Dems2)]
  # Nah, use original coding.
  DT[, votes.with.Dems2 := votes.with.Dems]

  # Plot tabulation by condition (Senators, Standard, and Info).
  cat("Tabulation of number votes with Dems by condition for
  those who cast 12 votes:\n")
  print(round(tab <- DT[tot.votes == 12,prop.table(xtabs(weight~votes.with.Dems2+assignedCondition),2)],3))

  pdf(sprintf("figures/%s",second.fig.name),width=10.5,height=6)
  par.old <- par(mar=c(4.1,4.1,1.1,0.6))
  # Combined barplot of Dem+Rep Senators and respondents
  # by condition.
  DT[,short.pid := ifelse(!is.na(sen.pid),sen.pid,substr(pid3_DtoR,1,1))]
  DT[short.pid != "O",condition2 := sprintf("%s-%s",short.pid,gsub("Info","Party split",gsub("Standard","Control",gsub("Senator","Senators",assignedCondition))))]
  DT[,condition3 := factor(condition2,levels=c("R-Senators","R-Control","R-Party split","D-Senators","D-Control","D-Party split"),ordered=T)]
  DT[,vwd3 := factor(2*(floor(votes.with.Dems2/2)),labels=c("0-1","2-3","4-5","6-7","8-9","10-11","12"))]

  # Plot tabulation by condition (Senators, Standard, and Info).
  cat("Tabulation of number votes with Dems by condition for
  those who cast 12 votes:\n")
  print(round(tab <- DT[tot.votes == 12,prop.table(xtabs(weight~votes.with.Dems2+assignedCondition),2)],3))

  # Combined density plot of Dem+Rep Senators and respondents
  # by condition.
  tab <- DT[tot.votes == 12 & short.pid != "O",prop.table(xtabs(weight~vwd3+condition3),2)]
  cols <- c(rev(brewer.pal(9,"Reds")[c(4,7,9)]),brewer.pal(4,"Blues")[2:4])
  ltys <- c(3,1,4,3,1,4)
  plot(x=c(1,nrow(tab)),y=range(tab),ann=F,axes=F,type='n')
  .x <- seq_len(nrow(tab))
  axis(1,at=.x,labels=rownames(tab))
  axis(2,las=2)
  title(main="", ylab="Proportion", xlab="Number votes with Democratic position")
  for (j in 1:ncol(tab)) {
    lines(x=.x,y=tab[,j],col=cols[j],type="l",lwd=4,pch=19,lty=ltys[j])
  }
  legend("top",legend=colnames(tab),col=cols,lwd=2,lty=ltys)
  par(par.old)
  dev.off()
}

#
# Three versions of figures.
#

# Weighted and Screened
.doFigures(T,T,"TreatmentEffectsByBillPid3-Screened.csv","Figure02InfoVsStndSupport.pdf","PreppedIRTData-Screened.csv","Figure06VotesWithDemsByCondition-v2.pdf")

# Weighted and Unscreened
.doFigures(T,F,"TreatmentEffectsByBillPid3-Unscreened.csv","FigureA04InfoVsStndSupport.pdf","PreppedIRTData-Unscreened.csv","FigureA05VotesWithDemsByCondition-v2.pdf")

# Unweighted and screened
.doFigures(F,F,"TreatmentEffectsByBillPid3-Screened-Unweighted.csv","FigureA10InfoVsStndSupport.pdf","PreppedIRTData-Unscreened.csv","FigureA11VotesWithDemsByCondition-v2.pdf")


#
# %%%%%%%%%%%%%%%%%%%%%%%%
# Loess of policy confidence on policy importance
# by issue.
# %%%%%%%%%%%%%%%%%%%%%%%%
#

# Data from AnalyzeSSI.do.
pol <- fread("RawData/ConfByPolicyImportance-Screened.csv")
# Map of confidence to importance to policy area.
map <- data.table(conf=c('conf_Q20_1',	'conf_Q20_2',	'conf_Q20_3',	'conf_Q20_5',	'conf_Q21_1',	'conf_Q21_2',	'conf_Q21_3',	'conf_Q21_4',	'conf_Q22_1',	'conf_Q22_2',	'conf_Q22_3',	'conf_Q22_5'), impt=c('impt_Q23_1',	'impt_Q23_2',	'impt_Q23_3',	'impt_Q23_5',	'impt_Q24_1',	'impt_Q24_2',	'impt_Q24_3',	'impt_Q24_4',	'impt_Q25_1',	'impt_Q25_2',	'impt_Q25_3',	'impt_Q25_5'), lab=c('Reducing gun violence',	'Making college affordable',	'Building and maintaining transportation system',	'Balancing the federal budget',	'Getting the unemployed back to work',	'Preventing gender discrimination in employment',	'Supporting domestic oil and gas production',	'Guaranteeing affordable housing',	'Guaranteeing access to affordable medical care',	'Making Medicare a program that works for patients',	'Promoting American export-oriented businesses',	'Reducing the number of abortions'))

.makeLines <- function(x,y,w=1,new.x) {
  # Make the predicted values of a smooth of y on x.
  # Arguments.
  #  x,y -- data points.
  #  w -- optional weights vector.
  #  new.x -- values of x at which to make predictions.
  l <- loess(y~x,weigths=w,span=2/3)
  return(predict(l,data.frame(x=new.x)))
}

# Data.table to hold predicted values.
new.x <- seq(0,3,.1)
res <- data.table(new.x=new.x)
# Loop over policy areas and make smooths.
for (i in 1:nrow(map)) {
  pred <- .makeLines(y=pol[[map[i,conf]]],x=pol[[map[i,impt]]],w=pol[["weight"]],new.x=new.x)
  # Only use values inside range of x.
  pred[new.x > max(pol[[map[i,impt]]],na.rm=T) | new.x < min(pol[[map[i,impt]]],na.rm=T)] <- NA
  res[,sprintf("pred%02d",i) := pred]
}

pdf("figures/FigureA07ConfVsImpLoess.pdf",width=10,height=6)
par.old <- par(mar=c(4.1,4.1,1.1,0.6),cex.lab=1.2,cex.axis=1.2)
# Matplot the lines.
matplot(x=res[,new.x],y=res[,sprintf("pred%02d",1:nrow(map)),with=F],type='l',lwd=3,col='black',lty=1,xlim=c(-0.05,3),ylim=c(-0.1,3),ann=F,axes=F)
axis(1,at=seq(0,3));axis(2,las=2,at=seq(0,3))
title(xlab="Importance of policy",ylab="Confidence distinguishing good from bad policies")

# Add means of all variables as cross hatches.
library(matrixStats)
abline(h=colWeightedMeans(x=as.matrix(pol[,map[,conf],with=F]),w=pol[,weight],na.rm=T),col='gray',lwd=2)
abline(v=colWeightedMeans(x=as.matrix(pol[,map[,impt],with=F]),w=pol[,weight],na.rm=T),col='gray',lwd=2)
# Add loess smooths back over top.
matlines(x=res[,new.x],y=res[,sprintf("pred%02d",1:nrow(map)),with=F],type='l',lwd=3,col='black',lty=1)
# Add text indicating distribution of responses.
tab <- prop.table(table(unlist(pol[,map[,conf],with=F])))
text(y=as.numeric(names(tab)),x=0,lab=sprintf("%1.0f%%",100*tab),pos=2)
tab <- prop.table(table(unlist(pol[,map[,impt],with=F])))
text(x=as.numeric(names(tab)),y=0,lab=sprintf("%1.0f%%",100*tab),pos=1)

par(par.old)
dev.off()
